function res = imageapprox(A, varargin)
% imageapprox     Make a (sparse) approximation of the image.
%                 The reconstructed image is: Ar = imagerestore(res);
% The main steps in this function are
% 1. The image is arraged into column vectors, X = myim2col(A, ...)
% 2. A sparse approximation of the columns in X is made using the
%    supplied dictionary or thresholding (only transform)
%    Sparseness is made based on targetPSNR or tsf. 
% 3. The weights 'Zw', and the DC component 'Zdc', may be quantized using a
%    uniform quantizer with step delta (and a threshold), uniquant.
% 4. 'Zw' and 'Zdc' may be rearranged into sequences and bit rate is then
%    estimated, or a byte sequence 'y' may even be made.
%
% The function returns the results in a struct. The input arguments are the
% image A and several pairs of  input arguments collected in a struct, a
% cell or as pairs: argName, argVal, ....  
%
% res = imageapprox(A, argName, argVal, ...);
%----------------------------------------------------------------------
%  res   a struct with the results, fields may be:
%        tr (transform), PSNR, PSNRbq (before quantizing), delta,
%        threshold, nonzeros, sf (sparseness factor), Ar (restored image),
%        Zdc (quantized DC component), Zw (sparse matrix with quantized
%        weights), Xdc (DC component before quantizing), W (sparse matrix
%        with the weights)
%
%  A     an image, a matrix of doubles (centered around 0), size Ma x Na.
%        If image size does not match transform size (dictionary) image
%        will be extended (by myim2col). A should be only one layer.
%        ex: A = double(imread('lena.bmp')) - 128; % lena.bmp is grayscale
%
%        The input arguments (options) may be:
%  'peak' the peak for image A (used in PSNR), default 255
%  't' or 'transform' is as in myim2col, default is 'none'. 
%  's' or 'size' is for the size of the transform, also as in myim2col.
%        Transform and its size may be overruled by the dictionary
%  'd' or 'dictionary' is the dictionary to use. It may be:
%        a. A dictionary do not need to be given, option may be omitted
%        b. The dictionary as a matrix D 
%        c. A struct with a field named 'D', and possible also 'transform'
%        d. The name of a mat-file, containing 'D' (and 'transform')
%        Cases c and d usually override 'transform' and 'size'
%
%        The sparse representation is done by sparseapprox.m, options used
%  'tPSNR' or 'targetPSNR' is target PSNR, default 0 (not used)
%        if tPSNR is not given 'tsf' may be used
%  'tsf' or 'targetSparsenessFactor' is target sparseness factor, it
%        includes DC componet (which is always selected). So when 8x8
%        blocks are used to form the columns in X, 1/64 < tsf < 1
%        
%        Quantizing is done by uniquant.m, options used
%  'del' is binsize in uniform quantizer. Use 0 for no quantizing, and
%        a negative value, ex -10, for adapting (using -del as first)
%  'thr' gives size of zero-bin from -thr to +thr, default is 'del'
%        which makes zero-bin twice the size of other bins.
%  'delDC' as 'del' but used for DC component, default 'del'
%  'thrDC' as 'thr' but used for DC component, defualt 'delDC'/2
%  'qLimit' gives limits for how much we want quantizing to reduce PSNR
%        given as [min, target, max], quantizing loop stops when PSNR is
%        reduced by d, and d is between min and max. May give
%        [min, max] and target is (min+max)/2, or just may give just
%        target, and min/max = target +- 0.05, default qLimit = 0.4
%  'eb' or 'estimateBits' estimates (or often finds) the number of bits
%        needed to code the quantized values, i.e. Zdc and Zw.
%        Value may be '' (default), 'none','Huff06', 'Arith06' or 'Arith07'
%        where empty string ('') not estimate bits and 'none' estimate bits
%        based on zero-order entropy but do no entropy coding
%
%  'v' or 'verbose' to indicate verboseness, default 0
%----------------------------------------------------------------------
% example:
%  A = double(imread('lena.bmp')) - 128;
%  B = double(imread('d:/bilder/USC_tiff/barbara.bmp')) - 128;
%  p1 = struct('dictionary','ex312Jan191930.mat', 'verbose',1); 
%  % without quantizing
%  r1 = imageapprox(A, p1, 'tPSNR',37.22);
%  r1 = imageapprox(A, p1, 'tsf',10192/numel(A));
%  r1 = imageapprox(A, 't','m79', 's',8, 'tPSNR',37, 'v',1);
%  r1 = imageapprox(A, 't','m79', 's',8, 'tsf',0.05, 'v',1);
%  % with quantizing
%  r1 = imageapprox(A, p1, 'tPSNR',37.6, 'del',18);
%  r1 = imageapprox(A, 't','m79', 's',8, 'tPSNR',37, 'del',-18, 'v',1);
%  p = struct('d','ex312Jan191930.mat', 'tPSNR',37, 'del',-18, 'v',1);
%  r1 = imageapprox(A, p, 'qLimit',[0.49,0.51], 'eb','Huff06');

%----------------------------------------------------------------------
% Copyright (c) 2009.  Karl Skretting.  All rights reserved.
% University of Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  04.12.2009  Made function
% Ver. 1.1  14.12.2009  Function almost ok. 
% Ver. 1.2  07.01.2010  Function made simpler, and it works fine
% Ver. 1.3  11.02.2010  Added som options
% Ver. 1.4  01.03.2010  Tested (but not throughoutly) and works ok.
%----------------------------------------------------------------------
% need functions: myimadjust, myim2col, mycol2im, myreshape (from DLE)
% and sparseapprox (from ICT)
% Vector Selection: the 'mpv2' Java package, see javaAccess.m
% from my Comp toolbox: uniquant, myreshape, mypred 

%% default options
peak = 255; 
tr = '';            % transform
Ms=0; Ns=0;         % size (of transform)
amet = 'extend';    % adjust method in myim2col (as in myimadjust)
tpsnr = 0;          % target PSNR
tsf = 0.1;          % or target sparseness factor
del = 0;
thr = 0;
deldc = 0;
thrdc = 0;
verbose = 0;
qlimit = [0.35, 0.4, 0.45];
eb = '';

%% get the options
nofOptions = nargin-1;
optionNumber = 1;
fieldNumber = 1;
while (optionNumber <= nofOptions)
    if isstruct(varargin{optionNumber})
        sOptions = varargin{optionNumber}; 
        sNames = fieldnames(sOptions);
        opName = sNames{fieldNumber};
        opVal = sOptions.(opName);
        % next option is next field or next (pair of) arguments
        fieldNumber = fieldNumber + 1;  % next field
        if (fieldNumber > numel(sNames)) 
            fieldNumber = 1;
            optionNumber = optionNumber + 1;  % next pair of options
        end
    elseif iscell(varargin{optionNumber})
        sOptions = varargin{optionNumber}; 
        opName = sOptions{fieldNumber};
        opVal = sOptions{fieldNumber+1};
        % next option is next pair in cell or next (pair of) arguments
        fieldNumber = fieldNumber + 2;  % next pair in cell
        if (fieldNumber > numel(sOptions)) 
            fieldNumber = 1;
            optionNumber = optionNumber + 1;  % next pair of options
        end
    else
        opName = varargin{optionNumber};
        opVal = varargin{optionNumber+1};
        optionNumber = optionNumber + 2;  % next pair of options
    end
    % interpret opName and opVal
    if (strcmpi(opName,'peak') && isnumeric(opVal))
        peak = opVal(1);
    end
    if (strcmpi(opName,'transform') || strcmpi(opName,'t'))
        tr = opVal;
    end
    if (strcmpi(opName,'size') || strcmpi(opName,'s'))
        if ((numel(opVal) == 1) && isnumeric(opVal))
            Ms = max(0, floor(opVal));
            Ns = Ms;
        elseif ((numel(opVal) == 2) && isnumeric(opVal))
            Ms = max(0, floor(opVal(1)));
            Ns = max(0, floor(opVal(2)));
        else
            error('imageapprox: illegal option size. We ignore it.');
        end
    end
    if (strcmpi(opName,'dictionary') || strcmpi(opName,'d'))
        if isnumeric(opVal)
            Ds = struct('D',opVal,'N',size(opVal,1),'K',size(opVal,2));
            % Ds.transform is not given (use tr)
        end
        if isstruct(opVal)
            Ds = opVal;
        end
        if ischar(opVal) 
            if (exist(opVal, 'file') == 2)
                Ds = load(opVal);
            else
                error(['imageapprox: can not read mat-file ',opVal,'.']);
            end
            if ~isfield(Ds,'D')
                error(['imageapprox: mat-file ',opVal,' has no dictionary D.']);
            end    
        end
        % if Ds exist it is a struct and contain at least field D
    end
    %
    if ( (strcmpi(opName,'targetPSNR') || ...
            strcmpi(opName,'tPSNR')) && isnumeric(opVal) )
        tpsnr = opVal(1);
    end
    if ( (strcmpi(opName,'targetSparsenessFactor') || ...
            strcmpi(opName,'tsf')) && isnumeric(opVal) )
        tsf = opVal(1);
    end
    if ( strcmpi(opName,'del') && isnumeric(opVal) )
        del = opVal(1);
    end
    if ( strcmpi(opName,'thr') && isnumeric(opVal) )
        thr = opVal(1);
    end
    if ( strcmpi(opName,'delDC') && isnumeric(opVal) )
        deldc = opVal(1);
    end
    if ( strcmpi(opName,'thrDC') && isnumeric(opVal) )
        thrdc = opVal(1);
    end
    if ( strcmpi(opName,'qLimit') && isnumeric(opVal) )
        if numel(opVal) == 1
            qlimit = [opVal-0.05, opVal, opVal+0.05];
        elseif numel(opVal) == 2
            qlimit = [opVal(1), (opVal(1)+opVal(2))/2, opVal(2)];
        elseif numel(opVal) > 2
            qlimit = [opVal(1), opVal(2), opVal(3)];
        end
    end
    if (strcmpi(opName,'estimateBits') || strcmpi(opName,'eb')) 
        if ischar(opVal)
            eb = opVal;
        else
            disp('imageapprox: illegal option estimateBits. We ignore it.');
        end
    end
    %
    if strcmpi(opName,'verbose') || strcmpi(opName,'v')
        if (islogical(opVal) && opVal); verbose = 1; end;
        if isnumeric(opVal); verbose = opVal(1); end;
    end
end

%% check and display arguments and options
if (size(A,3) > 1)
    error('imageapprox: not ready for several layers in image A.');
end
if ~(isa(A,'double'))
    error('imageapprox: image A is not double.');
end

if exist('Ds','var')   % else: do not use any dictionary
    if ~isfield(Ds,'D')   % Ds is a struct
        error('imageapprox: Missing dictionary, (no Ds.D)');
    end
    if isfield(Ds,'transform')
        % the former lifting schemme 'ks2' is replaced by 'm79'
        if strcmpi(Ds.transform,'ks2'); Ds.transform = 'm79'; end;
        %
        if (verbose && (length(tr) > 1) && ~strcmpi(Ds.transform,tr))
            disp(['Overrule transform: change from ',tr,' to ',...
                Ds.transform,'.']);
        end
        tr = Ds.transform;   % any given transform is overruled!
    end
    if (size(Ds.D,1) < 1)
        error(['imageapprox ERROR: size of dictionary is only ',...
            int2str(size(Ds.D,1)),'x',int2str(size(Ds.D,2)),'.']);
    end
    if (size(Ds.D,1) < 8)
        disp(['imageapprox WARNING: size of dictionary is only ',...
            int2str(size(Ds.D,1)),'x',int2str(size(Ds.D,2)),'.']);
    end
    if (size(Ds.D,1) == 64)
        Ms = 8; Ns = 8;
    end
end

if ((Ms == 0) || (Ns == 0))
    if strcmpi(tr,'none')
        Ms = 4; Ns = 4;
    else
        Ms = 8; Ns = 8;
    end
end

if ~(iscell(tr) || ischar(tr))
    error('imageapprox: transform (tr) is not char or cell.');
end
if strcmpi(tr,'');  % still default
    error('imageapprox: transform (tr) was not given.');
end

if (strcmpi(tr,'lot') || strcmpi(tr,'elt'))
    Ms = Ms - mod(Ms,2);
    Ns = Ns - mod(Ns,2);
end
if ( iscell(tr) || ...  % a wavelet lifting scheme may be given by a cell
        ( ischar(tr) && ... 
          ~strcmpi(tr,'none') && ...
          ~strcmpi(tr,'lot') && ...
          ~strcmpi(tr,'elt') && ...
          ~strcmpi(tr,'dct') )  )
    % a wavelet as in mylwt2
    Ms = max(2,2^floor(log2(Ms)));
    Ns = Ms;
end
%
if (del ~= 0)
    if (thr == 0); thr = del; end;
    if (deldc == 0); deldc = del; end;
    if (thrdc == 0); thrdc = deldc/2; end;
end
thr = abs(thr); 
deldc = abs(deldc);
thrdc = abs(thrdc);
% 
if verbose
    disp(' ');
    disp(['imageapprox: Input image is ',int2str(size(A,1)),...
        'x',int2str(size(A,2)),' of class ',class(A)]);
    disp(['Arguments are: transform = ',tr,...
          ', tPSNR = ',num2str(tpsnr),...
          ', peak (for PSNR) = ',num2str(peak),...
          ', tsf = ',num2str(tsf),...
          ', del = ',num2str(del),...
          ', thr = ',num2str(thr),...
          ', delDC = ',num2str(deldc),...
          ', thrDC = ',num2str(thrdc)]);
end
        
%% do the transform and make the column vectors from the image
[Ma, Na] = size(A);    % original size of image
X = myim2col(A, 't', tr, 's', [Ms, Ns], 'a',amet, 'n', [Ms, Ns], 'i', [Ms, Ns]);
[N,L] = size(X);
if (Ma*Na) == (N*L)
    imageadjusted = false;
    Maa = Ma;
    Naa = Na;
else
    imageadjusted = true;
    Maa = ceil(Ma/Ms)*Ms; % adjusted image size
    Naa = ceil(Na/Ns)*Ns;
    if (Maa*Naa) ~= (N*L)
        error('imageapprox: Can not match size of A and X.');
    end
end
if verbose
    disp(['Image A of size ',int2str(Ma),'x',int2str(Na),...
          ' (',int2str(Ma*Na),' pixels)', ...
          ' was divided into blocks of size ',int2str(Ms),...
          'x',int2str(Ns),' (transform: ',tr,'),']);
    disp(['and then arranged into X of size ',int2str(N),...
          'x',int2str(L),' (',int2str(N*L),' elements).']);
end

%%  Extract DC both for transform only and (transform and) dictionary
if strcmpi(tr,'none')  % extract DC as first row (or mean)
    Xdc = sqrt(N)*mean(X);
    Xr = ones(N,1)*(Xdc/sqrt(N));  % reconstructed using only Xdc
    X = X - Xr;
else
    Xdc = X(1,:);
    X(1,:) = zeros(1,L);
    Xr = [Xdc; zeros(N-1,L)];
end
%
Ar = mycol2im(Xr, 't', tr, 'imsize', [Maa, Naa], 's', [Ms, Ns]);
if imageadjusted
    Ar = Ar(1:Ma, 1:Na);
end
R = A - Ar;
PSNRdc = 10*log10((((Ma*Na)*peak^2)/sum(sum(R.*R)))); % using only DC
if (verbose) % display PSNR for image based on Xdc only
    disp(['Using only the ',num2str(L),...
        ' elements in the DC component, PSNRdc = ',num2str(PSNRdc),'.']);
end

%%  do a sparse representation of X using the dictionary or thresholding 
if (exist('Ds','var') == 1)    % sparse W, Xa = Ds.D * W,
    javaAccess();
    if (tpsnr > 0)
        % Kanskje har tSSE blitt fjernet fra sparseapprox JA
        tSSE = ((Ma*Na)*peak^2)*10^(-tpsnr/10);   
        if verbose
            disp(['imageapprox calls sparseapprox using', ...
                ' tSSE = ',num2str(tSSE),...
                ' (target PSNR = ',num2str(tpsnr),')']);
        end
        W = sparseapprox(X, Ds.D, 'javaORMP', 'tnz',N/2, 'tSSE',tSSE, ...
            'v',max(0,verbose-1));
        W = sparse(W);
    elseif (tsf < 1) %  
        tnnz = floor(Ma*Na*tsf) - L;
        if verbose
            disp(['imageapprox call sparseapprox using', ...
                ' tnz = ',int2str(tnnz),...
                ' (target sparseness factor = ',num2str(tsf),')']);
        end
        W = sparseapprox(X, Ds.D, 'GMP', 'tnz',tnnz, 'v',max(0,verbose-1));
        W = sparseapprox(X, Ds.D, 'javaORMP', 'tnz',sum(W~=0), ...
            'globalReDist',1, 'v',max(0,verbose-1));
        W = sparse(W);
    else 
        disp('imageapprox do not call sparseapprox (use W = D\X;).');
        W = Ds.D\X;    % method of frames (is NOT sparse)
    end
    S = sum(W ~= 0);  % selected number of non-zeros for each column
    sumS = sum(S);
    % find PSNRbq  (before quantizing)
    Xa = Ds.D*W;         % the sparse approximation (excluding DC)
    Ar = mycol2im(Xr+Xa, 't', tr, 'imsize', [Maa, Naa], 's', [Ms, Ns]);
    if imageadjusted
        Ar = Ar(1:Ma, 1:Na);
    end
    R = A - Ar;
    PSNRbq = 10*log10((((Ma*Na)*peak^2)/sum(sum(R.*R)))); % sparse rep
    if verbose    % display results so far
        disp(['After dictionary D is used: non-zeros including DC is ',...
            int2str(sumS+L),' and PSNR = ',num2str(PSNRbq,6) ]);
    end
else  % no dictionary, set Xa = W to the largest entries of X
    [temp,I] = sort(abs(X(:)),'descend');
    W = zeros(size(X));
    if (tpsnr > 0)
        tSSE = ((Ma*Na)*peak^2)*10^(-tpsnr/10);
        max_error = sum(sum(X.*X));  % = sum(temp.*temp)
        if (max_error > tSSE)
            temp(1) = max_error - temp(1).^2;  % sum of the rest squared
            if (temp(1) < tSSE);
                i = 1;
            else
                for i=2:numel(temp)
                    temp(i) = temp(i-1) - temp(i).^2;
                    if (temp(i) < tSSE); break; end;
                end
            end
        end
        W(I(1:i)) = X(I(1:i));  % use the largest
        if verbose
            disp(['Using only ',int2str(i),' largest of X in W,', ...
                ' tSSE = ',num2str(tSSE),...
                ' (target PSNR = ',num2str(tpsnr),')']);
        end
    elseif (tsf < 1) %
        tnnz = min(max(1, floor((Ma*Na)*tsf) - L), numel(X));
        W(I(1:tnnz)) = X(I(1:tnnz));  % use the largest
        if verbose
            disp(['Using only ',int2str(tnnz),' largest of X in W,', ...
                ' (target sparseness factor = ',num2str(tsf),')']);
        end
    else    
        W = X;
        if verbose
            disp('No sparse representation, (W = X).');
        end
    end
    Xa = W;
    S = sum(W ~= 0);  % selected number of non-zeros for each column
    sumS = sum(S);
    % find PSNRbq  (before quantizing)
    Ar = mycol2im(Xr+Xa, 't', tr, 'imsize', [Maa, Naa], 's', [Ms, Ns]);
    if imageadjusted
        Ar = Ar(1:Ma, 1:Na);
    end
    R = A - Ar;
    PSNRbq = 10*log10((((Ma*Na)*peak^2)/sum(sum(R.*R)))); % sparse rep
    if verbose    % display results so far
        disp(['After thresholding: non-zeros including DC is ',...
            int2str(sumS+L),' and PSNR = ',num2str(PSNRbq,6) ]);
    end
end

%% Quantizing and find restored image
adaptdelta = false;
if (del < 0)  % try to select an appropriate value
    del = -del;
    adaptdelta = true;
end
for dummycounter = 1:10 
    % in loop only display if (verbose >= 2)
    if (del > 0)
        Zdc = uniquant(Xdc, deldc, thrdc, 2000);
        Zw = uniquant(W, del, thr, 2000);
        Qdc = uniquant(Zdc, deldc, thrdc); % inverse quantizing
        Qw = uniquant(Zw, del, thr);   % inverse quantizing
        dcnz = full(sum(Zdc ~= 0));    % number of non-zeros in DC
        S = full(sum(Zw ~= 0));  % selected number of non-zeros for each column
        sumS = sum(S);
        if verbose >= 2
            disp(['Quantizing loop number ',int2str(dummycounter),':']);
            disp(['  Total number of non-zeros after quantizing ',...
                int2str(dcnz+sumS),' (Xdc:',int2str(dcnz),', W:',int2str(sumS),...
                '), sparseness ',num2str((dcnz+sumS)/numel(A)),'.']);
            disp(['  using del = ',num2str(del),...
                ', thr = ',num2str(thr),...
                ', deldc = ',num2str(deldc),...
                ', thrdc = ',num2str(thrdc),'.']);
        end
    else % no quantizing  (note Zdc and Zw are not defined)
        Qw = W;
        Qdc = Xdc;
        dcnz = sum(Xdc ~= 0);    % number of non-zeros in DC
        S = sum(W ~= 0);  % selected number of non-zeros for each column
        sumS = sum(S);
    end
    
    if strcmpi(tr,'none')  % find Xr from Qdc
        Xr = ones(N,1)*(Qdc/sqrt(N));  
    else
        Xr = [Qdc; zeros(N-1,L)];
    end
    if (exist('Ds','var') == 1)    % sparse W, Xa = Ds.D * W,
        Xa = Ds.D*Qw;       % the sparse approximation (excluding DC)
        Ar = mycol2im(Xr+Xa, 't', tr, 'imsize', [Maa, Naa], 's', [Ms, Ns]);
    else
        Ar = mycol2im(Xr+Qw, 't', tr, 'imsize', [Maa, Naa], 's', [Ms, Ns]);
    end
    if imageadjusted
        Ar = Ar(1:Ma, 1:Na);
    end
    R = A - Ar;
    PSNRq = 10*log10((((Ma*Na)*peak^2)/sum(sum(R.*R)))); % after quant.
    
    if verbose >= 2
        if (del > 0)
            disp(['  Quantizing reduced PSNR with ',...
                   num2str(PSNRbq-PSNRq),' dB to PSNRq = ',num2str(PSNRq),'.']);
        end
        disp(['  Sparseness factor is ',int2str((dcnz+sumS)),...
            '/',int2str(numel(A)),' = ',num2str((dcnz+sumS)/numel(A)),'.']);
    end
    if (~adaptdelta); break; end;
    ratios = [thr, deldc, thrdc, del*del]/del;
    d = (PSNRbq-PSNRq); % difference (=reduction) in PSNR due to quatizing
    if ((qlimit(1) < d) && (d < qlimit(3))); break; end;
    % what the factor f should be is here ad-hoc, this often work 
    % (adapted to m79 transform)
    f = interp1([0   0.2  0.3  0.4  0.6   3    10   50],...
                [2.5 1.5  0.5  0.3  0.1 0.08 0.09  0.1],d);
    del = 10^(f*(qlimit(2)-d)) * del;
    if verbose >= 2
        disp(['  del was ',num2str(ratios(4)),...
              ', d (PSNRbq-PSNRq) = ',num2str(d), ...
              ', f = ',num2str(f),' qlimit(2) = ',num2str(qlimit(2))]);
        disp(['  In the end of quantizing loop ',int2str(dummycounter),...
              ' del was changed to ',num2str(del)]);
    end
    thr = del*ratios(1);
    deldc = del*ratios(2);
    thrdc = del*ratios(3);
end
if (verbose == 1) && (del > 0)
    disp(['Quantizing using del = ',num2str(del),...
        ', thr = ',num2str(thr),...
        ', deldc = ',num2str(deldc),...
        ', thrdc = ',num2str(thrdc),'.']);
    disp(['  ==>  nonzeros = (',int2str(dcnz),'+',int2str(sumS),...
        ') = ',int2str(dcnz+sumS),...
        ', sparseness = ',num2str((dcnz+sumS)/numel(A)),...
        ' and PSNR = ',num2str(PSNRq),'.']);
end
if (verbose >= 1) && (del == 0)
    disp('No quantizing was done.');
end

%% reshape and entropy coding
bits = [];
ebr = [];
y = [];
if ((numel(eb) > 1) && exist('Zdc','var') && exist('Zw','var'))   
    % estimate bits
    if (exist('Ds','var') == 1)  % a dictionary was used
        xCw = myreshape(Zw, 'method',2, 'v',0);
        xCdc = mypred( reshape(Zdc,Maa/8,Naa/8), 'nofS',3, 'v',0 );
        xC = cell(numel(xCw)+numel(xCdc),1);
        % disp([numel(xCw), numel(xCdc)]);
        for i=1:numel(xCw); xC{i}=xCw{i}; end;
        for j=1:numel(xCdc); xC{i+j}=xCdc{j}; end;
        [bits, ebr, y] = estimateBits(xC, eb, verbose-1);
    else  % only transform
        if (sum(abs(Zw(1,:)))==0) && (size(Zw,2)==numel(Zdc))
            % put DC into W
            Zw(1,:) = Zdc;
            xC = myreshape(Zw, 'method',1, 'sortrows',1, 'v',0);
            Zw(1,:) = 0;  % and remove it again
        else
            xC = myreshape(Zw, 'method',1, 'sortrows',1, 'v',0);
            xC{end+1} = [Zdc(1), diff(Zdc)];
        end
        [bits, ebr, y] = estimateBits(xC, eb, verbose-1);
    end
    if verbose
        disp(['Bit estimate (by ',eb,') is ',int2str(bits),...
            ' giving PSNR = ',num2str(PSNRq,'%5.2f'),...
            ', bit rate =',num2str(bits/(Ma*Na)),'.']);
    end
end

if ~exist('Zdc','var'); Zdc = []; end;
if ~exist('Zw','var'); Zw = []; end;
if ~exist('Ds','var'); Ds = []; end;
res = struct( 'text', ['Results returned by imageapprox.m, ',datestr(now())], ...
              'tr', tr, ...
              'sizes', struct('Ma',Ma,'Na',Na,'Maa',Maa,'Naa',Naa,...
                              'Ms',Ms,'Ns',Ns,'N',N,'L',L), ...
              'peak',peak, ...                
              'PSNR', PSNRq, ...
              'PSNRbq', PSNRbq, ...
              'PSNRdc', PSNRdc, ...
              'del', del, ...
              'thr', thr, ...
              'deldc', deldc, ...
              'thrdc', thrdc, ...
              'qlimit', qlimit, ...
              'nonzeros', (dcnz+sumS), ...
              'sf', (dcnz+sumS)/numel(A), ...
              'Ar', Ar, ...
              'Zdc', Zdc, ...
              'Zw', sparse(Zw), ...
              'bits', bits, ...
              'ebr', ebr, ...
              'ecm', eb, ...
              'y', y, ....
              'Qdc', ' Restored DC component: Qdc = uniquant(Zdc, deldc, thrdc); ', ...
              'Qw', ' Restored weights: Qw = uniquant(Zw, del, thr); ', ...
              'W', sparse(W), ...   
              'Ds', Ds );  
                  
return

